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Synthetic  3-D  Atmospheric  Temperature  Structure:  A  Modei 
for  Known  Geophysical  Power  Spectra  Using  a  Hybrid 

Autoregression  and  Fourier  Technique 


1.  INTRODUCTION 


Atmospheric  fluctuations  in  wind  speed,  temperature,  and  density  are  characterized  by 
continuous  power  spectral  density  functions.  For  example,  one-dimensional  horizontal  wind 
speed  PSDs  are  found  to  have  horizontal  wavenumber  log-log  slopes  of  about  -5/3  and  vertical 
wavenumber  log-log  slopes  of  -3.  Such  spectra,  parameterized  by  spectral  slope,  correlation 
lengths,  and  moments  of  probability  density  functions,  are  often  used  in  simulating  an 
environment  or  predicting  atmospheric  structure.  Multidimensional  fast  Fourier  transform 
synthesis  provides  a  means  for  filtering  white  noise  with  spatial  filters  to  simulate  a 
stationary  time  or  spatial  data  set.  Autoregression  synthesis  provides  a  fast  means  for 
simulating  a  one-dimensional  non- stationary  spatial  structure  sequence. 

The  Phillips  Laboratory  Strategic  High  Altitude  Atmospheric  Radiance  Code  (SHARC)^  uses 
first  principles  to  calculate  point  to  space  and  limb  viewing  atmospheric  background  infrared 
(IR)  radiance  and  transmittance  under  both  local-thermal-equilibrium  (LTE)  and  non-local- 


(Received  for  publication  14  June  1994) 
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Radiance  Code,  Phillips  Laboratory  Tech.  Rpt.,  PL-TR-9 1-2071.  ADA239008 
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thermal-equilibrium  (NLTE)  conditions  above  50  km.  Release  3  of  the  SHARC  code^  predicts  IR 
radiation  and  transmittance  in  the  1-40  pm  spectral  region  and  includes  important  bands 
from  the  major  isotopes  of  NO,  CO,  H2O,  O3,  OH,  CO2.  CH4,  and  NO"*".  Specific  local  atmospheric 
environments  can  be  specified  through  region  definitions,  and  diurnal  characteristics  can  be 
specified  through  user  or  program  generated  multiple  vertical  concentration  profiles.  A 
subroutine  module  breaks  a  given  line-of-sight  (LOS)  specification  into  small  segments  and 
determines  the  composition  and  properties  of  each  segment.  Each  segment  is  determined  by 
the  intersection  of  the  LOS  with  an  altitude  layer  boundary,  defined  in  input  atmospheric 
profiles.  Appropriate  profiles  of  temperature,  pressure,  and  molecular  state  densities  are 
determined  for  each  segment. 

A  future  release  of  SHARC  will  have  the  ability  to  provide  realizations  of  atmospheric 
infrared  volume-emission  perturbations  that  occur  from  fluctuations  in  temperature  and 
density  of  the  contributing  molecular  species.  Version  4  of  the  SHARC  code  envisions  a 
capability  to  evaluate  radiance  structure  from  estimated  variances  in  the  standard 
temperature  and  density  profiles.  The  algorithms  will  simulate  IR  fluctuations  that  must 
depend  on  relatively  small  fluctuations  in  atmospheric  species  number  densities,  vibrational 
state  populations,  and  the  kinetic  temperatures  along  a  given  line-of-sight.^  "^  Where  NLTE 
effects  dominate,  (generally  above  50  km)  a  small  fluctuation  in  kinetic  temperature  can 
produce  correlated  changes,  anti-correlated  changes,  or  no  change  in  the  vibrational  state 
temperature.  Such  changes  ultimately  depend  on  the  relative  contributions  from  total  number 
density,  temperature-dependent  kinetic  rates,  and  radiative  relaxation.  A  proper  description  of 
the  temperature/density  field  is  thus  needed  to  enable  SHARC  to  correctly  compute  the 
radiance  structure  field. 

To  provide  a  realistic  but  practical  two-dimensional  structure  scene  capability  requires 
creative,  efficient,  and  tested  algorithms.  This  report  presents  a  method  of  producing 
stochastic  three  dimensional  non- stationary  synthetic  spatial  structure  from  a  hybrid  one¬ 
dimensional  (vertical)  autoregressive  component  and  a  two-dimensional  (horizontal)  Fourier 


2  Gruninger,  J.,  Sundberg,  R.L.,  Duff,  J.W.,  Bernstein,  L.S.,  Matthew,  M.W.,  Adler-Golden,  S., 
Robertson,  D.,  Sharma,  R.,  Brown,  J.H.,  Healey,  R.,  and  Vail  J.  (1994)  SHARC  -  3,  A  model  for 
infrared  radiance  at  high  altitudes,  Proc.  SPIE  -  The  International  Society  of  Optical 
Engineering,  V;2223,  Orlando,  Florida. 

3  Sundberg,  R.L..  Gruninger,  J.,  De,  P.,  Brown,  J.H.  (1994)  Infrared  radiance  fluctuations  in 
the  upper  atmosphere,  Proc.  SPIE  -  The  International  Society  of  Optical  Engineering,  V:2223, 
Orlando,  Florida. 

^  Sears,  R.D.,  Strugala,  L.A.,  Newt,  J.,  Robertson,  D.,  Brown,  J.H.,  and  Sharma,  R.  (1994) 
Simulation  of  the  infrared  structured  earthlimb  background  using  the  SHARC  radiance  code, 
32nd  Aerospace  Sciences  Meeting  and  Exhibit,  Reno,  NV. 
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transform  component.  In  this  report  we  extend  the  results  expressed  in  two  previous  reports 
that  dealt  with  one-dimensional  autoregression®  and  two-dimensional  autoregression/moving 
average®  structure  simulation. 


2.  GEOMETRY 


SHARC  is  capable  of  predicting  earthlimb  radiance  along  a  sensor  line-of-sight  (LOS)  as  a 
function  of  LOS  and  tangent  altitude.  In  this  report,  temperature  and  density  perturbations  are 
incorporated  on  the  mean  atmosphere  to  enable  SHARC  calculations  of  stochastic  radiance 
structure.  An  appropriately  configured  three-dimensional  database  of  stochastic  temperature 
structure  is  needed.  To  fill  this  need,  the  non-statlonary  database  has  been  developed  in  a 
three-dimensional  matrix  that  is  described  geometrically  by  a  segment  of  a  sphere.  The 
spherical  segment  consists  of  a  volume  70  km  in  radial  height  by  approximately  3000  km  long 
by  approximately  51.2  km  wide.  The  horizontal  matrix  elements  lie  in  a  spherical  surface  at 
radial  distance  “r”.  Since  we  wish  the  database  (that  is,  the  synthetic  atmosphere)  to  be 
applicable  for  observation  by  a  modeled  sensor,  we  orient  the  longer  horizontal  dimension 
along  the  sensor  LOS  direction  and  the  shorter  horizontal  dimension  along  the  direction 
transverse  to  the  line-of-sight.  The  spatial  resolution  of  a  matrix  element  has  an  incremental 
altitude  of  0. 1  km  by  an  incremental  transverse  distance  of  approximately  0. 1  km  by  an 
incremental  LOS  distance  of  approximately  15.7  km.  Thus,  the  database  matrix  consists  of 
701  vertical  elements  aligned  along  radial  lines,  512  transverse  horizontal  elements,  and  192 
line-of-sight  horizontal  elements.  The  actual  data  spacing  in  the  horizontal  dimension  must 
Increase  slightly  with  altitude  so  that  points  in  the  vertical  dimension  stay  on  radial  lines. 
The  altitude  range  accounts  for  dominant  structure  effects  (50-120  km).  Also  its  vertical  and 
transverse  dimensions  are  sufficient  to  fill  the  field  of  a  limb-viewing  two-dimensional  sensor 
array.  The  vertical  and  transverse  spatial  resolution  are  governed  by  the  vertical  and 
horizontal  correlation  lengths  of  the  temperature  variations  but  also  are  chosen  to  provide 
high  field-of-view  spatial  resolution.  The  geometrical  LOS  dimensions  are  chosen  to  provide 
sufficient  length  to  describe  the  integrated  path  radiance.  Figure  1  illustrates  the  spherical 
geometry  and  Table  1  summarizes  the  database  specifications. 


®  Brown,  J.H.  (1993)  Atmospheric  Structure  Simulation:  An  Autoregressive  Model  for  Smooth 
Geophysical  Power  Spectra  with  Known  Autocorrelation  Function,  Phillips  Laboratory  Tech. 
Rpt,  PL-TR-93-2185,  ERP  No.  1128,  ADA  276691. 

®  Brown,  J.H.  (1993)  Atmospheric  Structure  Simulation:  An  ARMA  Model  for  Smooth  Isotropic 
Two-Dimensional  Geophysical  Power  Spectra,  Phillips  Laboratory  Tech.  Rpt.,  PL-TR-93-2224, 
ERP  No.  1 132,  ADA  280476. 
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a  b 


Figure  1.  Database  Geometry.  Figure  la  illustrates  the  volume  of  the  spherical 
segment  defined  by  the  database.  The  radial  altitude  range  starts  at  50  km  and 
ends  at  120  km.  At  50  km  altitude  the  curved  length  in  the  LOS  direction  measures 
3014  km  while  the  curved  length  in  the  transverse  direction  measures  51.2  km. 
These  distances  increase  slightly  with  altitude  due  to  the  curved  geometry.  Figure  lb 
illustrates  the  resolution  volume  at  50  km.  The  radial  altitude  resolution  is  0. 1  km 
while  the  resolution  of  the  curved  lengths  in  the  LOS  and  transverse  directions  is 
15.7  km  and  0.1  km  respectively.  Again  the  horizontal  resolution  lengths  increase 
slightly  with  altitude. 


Table  1.  Database  Specifications 


Direction 

N  pts 

Ax 

(km) 

N  A  X 
(km) 

freq.  res. 

1/N  A  X  (cy/km) 

Nyquist  freq. 
1/2  A  X  (cy/km) 

Vertical 
(50-120  km) 

700 

0.1 

70 

Transverse 

512 

0.1 

51.2 

0.0195  (HIGH) 

5 

LOS 

192 

15.7 

3014. 

0.00033 

0.0318  (LOW) 
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3.  AUTO-REGRESSION  THEORY 


This  section  briefly  reviews  autoregressive  digital  spectral  estimation  as  presented  by  S. 
Lawrence  Marple^  and  Steven  M.  Kay.® 

An  autoregressive  moving  average  (ARMA)  model  for  a  discrete  spatial  series  vfnj,  that 
approximates  deterministic  and  stochastic  processes  can  be  represented  by  the  filter  linear 
difference  equation: 


v(n) 


-  ia(j)  v(n-j)+  ib(j)e(n-j) 

j=i  J-o 


in  which  vfrfl  is  the  output  sequence  and  t(n)  a  white  noise  input  driving  sequence.  The  a(}) 
and  blJ)  coefficients  form  the  autoregressive  and  moving  average  portions  of  the  ARMA  model 
respectively.  A  z-transform  analysis  of  the  difference  equation  shows  that  the  ARMA  power 
spectral  density  (Parma^^^^  i®' 


^ARMA 


(k)  =  Ava^ 


m 

A(k) 


for  -oo  <  fc  <  oo 


where  k  is  the  spatial  frequency  and  where. 


p 

A{k)  =  1+2  o.{J)  exp  {-2nikjAv) 
J=i 


and, 


B{k)  =  1  +  2  b(j)  exp  {-2nikJAv),  with  b(0)  =  1. 
j=i 


7  Marple,  S.L.  (1987)  Digital  Spectral  Analysis  with  Applications,  Chapter  6,  Prentice-Hall, 
Englewood,  New  Jersey. 

®  Kay,  Steven  M.  (1988)  Modern  Spectral  Estimation,  Theory  &  Application,  Prentice-Hall, 
Englewood,  New  Jersey. 
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Here  Av  is  the  sampling  interval,  cr^  is  the  variance  of  the  white  noise  process,  and  i  =  If 
all  the  moving  average  coefficients  are  zero,  except  b(0)  =  1,  then. 


v(n)  =  -  i  a{j)v{n  -  j)  +  e(n) 
J=i 


and  the  process  is  strictly  autoregressive  of  order  p.  The  autoregressive  PSD  becomes: 


2Ava^ 

|A(/c)P 


=  2Av  X  ryy(j)exp(-27ci/gAv),  for  0  < /c<  oo 
J=-“ 


where  is  the  real  and  even  autocorrelation  sequence  at  j,  and  is  the  variance  of  e(n).  The 
relationship  between  the  autocorrelation  sequence  and  the  pure  autoregressive  model  is: 


rw(m) 


p 

-  J) 

J=1 

p 

-'Z^uyyvi-j)  +  0% 

j=i 


for  m  >  0 
for  m  =  0 
for  m  <  0 


This  expression  may  be  evaluated  for  the  p  +  1  lag  indices  0  <m<p  by: 


'rw(O)  • 

•  r^(-p) ' 

•  r^(a(l)) 

(  ^  ' 

a(l) 

II 

f„2  ^ 

0 

f  wip) 

Mp). 

.0. 
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This  expression  forms  the  autoregressive  Yule-Walker  equations.  Given  the  autocorrelation 
sequence  for  lags  0  to  p,  the  autoregressive  coefficients  may  be  found  from  the  above.  Since 
r^[-k)  =  r*^(k),  the  autocorrelation  matrix  is  both  Toeplitz  and  Hermitian.  A  standard 
“Levinson”  algorithm  that  takes  advantage  of  the  Hermitian-Toeplitz  matrix  equation  was 
employed  to  solve  for  the  AR  parameters.  The  same  Yule-Walker  equations  also  occur  if  we 

attempt  to  solve  the  problem;  find  the  “best”,  in  a  least  square  sense,  set  of  equations  that 

P 

determine  the  coefficients  that  predict  v  from  v(n)  =  -  X  ci(j)v(n  -j).  where 

j=i 

(v(n)  -  v(n))^  =  .  In  summary,  the  Levinson  recursion  computes  sets  of  coefficients 


1*^1  (l)>  |>  1^2  (^)’  *^2 ‘  ’  l^p  (p)’  *^p  } 


where  the  final  set  at  order  p  is  the  desired  solution  of  the  Yule-Walker  expressions.  For  the 
AR(p)  process,  a^ij]  =  a[j]  forj  =  1,2,3,..,  p  and  is  the  minimum  prediction  error,  that  is, 

Op  =  o^  =a^n  =^[v*[n](v[n])-v[n]].  The  algorithm  is  initialized  by: 


rw[l] 

w[0] 


af  =  (l-|ai[lf]r^[0] 


with  the  recursion  for  J  =  2,3,..,  p  given  by; 


a/[j]  =  - 


j-i 

rw[j]+  X  aj_i[p]rvv[j-li] 
_ iizl _ 

‘’5-1 


“  J  h]  =  o.j-\  h]  +  a  J  [  j]  a*_i  [  J  -  Ti]  q  =  1, 2, . . . ,  J  - 1 

a5  =  (l-|aj  [j]\^)o%. 
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®  Tatarski,  V.I.  (1961)  Wave  Propagation  in  a  Turbulent  Medium,  McGraw-Hill. 

Futterman,  W.I.,  Schweitzer,  E.L..  and  Newt.  J.E.  (1991)  Estimation  of  scene  correlation 
lengths,  characterization,  propagation,  and  simulation  of  sources  and  backgrounds,  Proc. 
SPIE  -  The  International  Society  of  Optical  Engineering,  V:  1486,  pp  127- 140,  Orlando,  Florida. 

11  For  the  integral,  see  for  example,  Gradshte)m,  I.S.  and  Ryzhik,  l.M.  (1965)  Table  of 
Integrals  Series  and  Products,  eq  3.241.4,  Academic  Press. 
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The  relationship  between  the  frequency  domain  PSD  and  the  time  or  spatial  domain 
autocorrelation  function  is  specified  by  their  Fourier  transform  pairs.  Thus  the 
autocorrelation  function,  (ACF),  for  the  real  even  PSD  function  is: 


ACF{x)  =  IFFT  (PSD)  =  2 


I 


cos{2nxk)dk 

2) 


1 

u+— 


0  +  2 


The  integral  is  solved  in  terms  of  the  Bessel  function  of  the  second  kind  of  fractional  order, 


as: 


•  12 


ACF(x)  = 


g2  ^(1-d)  (27cax:)^  (2Tcax) 

r(u) 


The  parameter  "a”  can  be  expressed  in  terms  of  the  integral  scale  or  equivalent  width,  L,  of 
the  autocorrelation  function.  An  equivalent  width  is  defined  as  the  area  of  a  function  divided 
by  its  central  ordinate,  or. 


{  fix) dx 
L  =  ^ - 

/(O) 


Thus  the  large  scale  correlation  length  is  defined  by  integrating  ACF  (x)  over  all  positive 
values  of  x. 


[“ACF(x)dx 

=:i0 - - - 


12  Gradshteyn,  I.S.  and  Ryzhik,  l.M.  (1965)  Table  of  Integrals  Series  and  Products,  eq  8.432.5, 
Academic  Press. 

13  Bracewell,  R.N.  (1978)  The  Fourier  Transform  and  Its  Applications,  McGraw-Hill,  New  York. 
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where  is  the  autocorrelation  at  zero  separation.  Since 

f  ACF  {x)dx  =  —  f  ACF  {x)dx  —  =  PSD{0)  and  ACF  (0)  =  the  equivalent  width  or  cor- 

JO  O  J  —  oo  O 


JO  2  2 

relation  length  as  a  function  of  the  PSD  may  also  be  written  as: 


PSD  jO)/ 2  _  PSD(0) 

PSD{k)dk 

J  — oo 


so  that, 


a^r(u+-]a^‘'  r(o+- 

—  _ V  ^  ^  ^  ^  QJ- 

^  2a^  2-\fKr{v)a 


1  2V^r(u) 

a  -  ,  V  . 

r  v+— 

<  2j 


5.  SIMULATION  OF  CORRELATED  VERTICAL  COMPONENTS 


In  practice,  a  given  one-dimensional  single-sided  PSD  was  constructed  from  the  PSD  and 
ACF  power  law  models: 


PSD(fc)  = 


2afa^'’rlu+- 

^  2 

l  2j 

20;^  Av 

V7cr(u)|a^-i-lc^|  2 

N 

,  V  i2nkAv 

l  +  2^aje 

>1 

,  for  0  <  k<°° 


ACFIv)  =  {2nayy  K^{2nav) 

r{v) 
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where, 


=  residual  error 
V  =  vertical  distance 
Av  =  vertical  spacing 
Uj  =  model  parameters 

i  = 

/c=  frequency  (km'l) 


Conceptually,  radial  line  segments  having  an  origin  at  the  center  of  the  earth  define  a 
segment  of  a  sphere.  Along  the  geometrical  set  of  these  radial  or  vertical  line  segments,  a 
spatial  data  sequence  was  simulated  using  a  normally  distributed  set  of  pseudorandom 
numbers  G(J).  In  this  case,  there  were  512x192  independent  sets  of  random  numbers 
corresponding  to  vertical  elements  along  the  radial  lines.  The  sequence,  G(  J),  was  initialized 
by  an  arbitrary  seed  and  thus  the  set  of  random  numbers  differed  for  each  line.  G(  J)  consisted 
of  a  mean  equal  to  zero  and  variance  equal  to  1.  The  simulation  of  701  sequence  values  along 

each  vertical  (radial)  line  was  generated  from  six  autoregression  coefficients  by  the  expression: 
6 

Y{J)  =  G{J)-'^a  where  the  values  were  determined  from  the  previously  described 

J=i 

Levinson  algorithm.  Note  that  the  first  level  of  non-stationarity  is  built  into  the  database  by 
allowing  the  aj  values  to  dynamically  change  with  altitude  as  a  function  of  the  vertical 

correlation  length.  The  first  51  values  of  Y(J)  for  each  radial  line  were  discarded  to  allow  for 
filter  relaxation.  In  practice,  for  altitudes  above  50  km,  the  calculational  scheme  stepped 
along  at  each  altitude  layer.  First  the  vertically  correlated  component  for  a  particular  layer 

o 

(the  top  layer  in  the  array)  was  computed  with  a  =1.  The  layer  was  then  copied  to  a  second 
array  where  the  appropriate  horizontal  correlation  and  o^  were  applied  for  that  layer.  Since 
and  horizontal  correlation  length  change  with  altitude,  these  parameters  were  applied  to  the 
horizontal  layer  in  the  copied  cirray  to  account  for  the  second  level  of  non-stationarity.  Thus 
each  layer  in  the  copied  array  contained  both  the  vertical  and  horizontal  correlations  before 
the  computations  proceeded  to  the  next  higher  layer.  The  scheme  is  represented  in  Figure  2. 
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The  altitude  dependence  of  the  relative  temperature  variance,  horizontal  correlation 
length,  and  vertical  correlation  length  are  adopted  from  Strugala,  et.al. Their  values  and 
our  fitted  models  of  them  are  shown  in  Figures  3-5.  Logistic  dose  response  equations  were  used 
in  fitting  the  data.  The  outside  curves  in  each  graph  represent  the  95  percent  confidence  limits. 
The  fitted  curves  were  used  to  evaluate  the  parameters  as  a  function  of  altitude. 


Figure  2.  Computational  Scheme.  First,  correlations  along  the  vertical  lines  are 
calculated  (one  layer  at  a  time)  and  then  the  horizontal  correlation  calculation 
is  applied  to  each  layer. 


14  Strugala,  L.A.  (1991)  Development  of  high  resolution  statistically  non-stationaiy  infrared 
earthlimb  radiance  scenes,  characterization,  propagation,  and  simulation  of  sources  and 
backgrounds,  Proc.  SPIE  -  The  International  Society  for  Optical  Engineering,  V:1486,  pp.  176-187, 
Orlando,  Florida. 

15  Strugala,  L.A.  (1993)  Production  of  statistically  nonstationaiy  stochastic  structure 
realizations  for  infrared  background  scene  simulations.  Optical  Engineering,  V32(No.  5)., 
pp.  993-1001. 
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Horizontal  Correlation  Length  (km) 

Eqn  8013  y=a+b/(1 +(x/c)d)  [LogisticDoseRsp] 
r2=0.991 459845  DF  Adj  r2=0.9891 82471  FitStdErr=2.07219665  Fstat=61 9.1 67443 

a=32.764354  b=52.968001 
C=1 35.0791 3  d=-6.4666086 


Rank  18  Eqn  8013  y=a+b/(1+(x/c)d)  [LogisticDoseRsp] 


Coef  Det 
0.9914598452 


DF  Adj  r2 
0.9891824705 


Fit  Std  Err 
2.0721966466 


F-value 

619.16744310 


Farm  Value 
a  32.76435426 

b  52.96800100 

c  135.0791315 

d  -6.46660858 


Std  Error 
1.295628218 
2.144297606 
2.234742679 
0.680074818 


t-value 

25.28839199 

24.70179552 

60.44504933 

-9.50867230 


95%  Confidence  Limits 


30.01210311 

48.41295545 

130.3319571 

-7.91126428 


35.51660541 

57.52304655 

139.8263059 

-5.02195287 


Figure  3.  Horizontal  Correlation  Lengths  Plotted  as  a  Function  of  Altitude. 
A  logistic  dose  response  equation  is  fitted  to  the  data.  The  outside  curves 
represent  the  95  percent  confidence  limits. 
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Vertical  Correlation  Length  (km) 

Rank  4  Eqn8013  y=a+b/(1+(x/c)^)  [LogisticDoseRsp] 
r2=0.977486661  DF  Adj  r2=0.9721 89405  FitStdErr=0.546523677  Fstat=260.508669 

3=1.7292654  6=8.3478628 
c=1 61. 52008  d=-7.9024983 


Rank  4  Eqn8013  y=a+b/(1+(x/c)'^)  [LogisticDoseRsp] 


Coef  Det 
0.9774866610 


DF  Adj  r2 
0.9721894048 


Fit  Std  Err 
0.5465236769 


F-value 

260.50866885 


Farm  Value 
a  1.729265412 

b  8.347862811 

c  161.5200848 

d  -7.90249825 


Std  Error 
0.210585414 
0.536310568 
3.853308016 
1.339177504 


t-value 

8.211705540 

15.56535207 

41.91725244 

-5.90100881 


95%  Confidence  Limits 


1 .286422835 
7.220048927 
153.4169175 
-10.7186706 


2.172107988 

9.475676694 

169.6232522 

-5.08632592 


Figure  4.  Vertical  Correlation  Lengths  Plotted  as  a  Function  of  Altitude.  A 
logistic  dose  response  equation  is  fitted  to  the  data.  The  outside  curves  represent 
the  95  percent  confidence  limits. 
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6T/T 


(5T/T)rms  Values 
[LgstcDoseRsp]  y=a+b/(1+(x/c)^) 

r2=0.99459057  DF  Adj  r2=0.99026303  FitStdErr=0.0097220158  Fstat=367.72483 
3=0.0087695291  b=0.31 03953 
c=1 18.04806  d=-3. 1509968 


Rank  931  Eqn8013  [LgstcDoseRsp]  y=a+b/(1+(x/c)'^) 


Coef  Det 
0.9945905715 


DF  Adj  r2 
0.9902630287 


Fit  Std  Err 
0.0097220158 


F-value 

367.72482563 


Farm  Value 
a  0.008769529 

b  0.310395297 

c  118.0480550 

d  -3.15099683 


Std  Error 
0.009675816 
0.024196036 
5.745293881 
0.473460616 


t-value 

0.906334825 

12.82835344 

20.54691326 

-6.65524591 


95%  Confidence  Limits 


-0.01500493 

0.250943187 

103.9312859 

-4.31433753 


0.032543990 

0.369847408 

132.1648241 

-1.98765613 


Date  Time  File  Source 

Jun  1,1994  10:59:59  AM  d:\tcwin2\rms1.prn 


Figure  5.  Root  Mean  Square  Values  of  the  Relative  Temperature  Fluctuations 
Plotted  as  a  Function  of  Altitude.  A  logistic  dose  response  equation  is  fitted  to 
the  data.  The  outside  curves  represent  the  95  percent  confidence  limits. 
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After  simulating  the  vertical  spatial  sequences  of  a  given  PSD,  the  PSD  of  the  simulated 
vertical  series  was  checked  using  a  forward  and  backward  estimation  method.  A  set  of  6  bj 
coefficients  were  found  to  estimate  the  simulated  PSD  from  the  formula, 


PSD^,,{k)  = 


2o^  Av 


1  + 

J=i 


i2nkAv 


2  ’ 


where  is  described  below  and  where  Av  is  the  same  spacing  used  to  generate  the  synthetic 
series  (0. 1  km). 

For  a  given  altitude,  the  forward  and  backward  method  solves  the  following  least  squares 

problem.  Given  13  discrete  vertical  series  values  Y(J),  for  J  =  1,  2 .  13,  and  where  J=7  is  the 

layer  in  question,  we  find  bj  values  that  minimize  ERR,  where. 


( 

13 

f  6 

2^ 

( 

7 

^6  ^ 

2" 

ERR  = 

z 

Y(J)-^bjY(J-j) 

+ 

z 

Y(J)-£bj  Y(J+j) 

J=7 

<  >1  > 

J 

J=1 

V 

.  >1  J 

y 

ERR 

Then  is  equal  to  - .  In  practice,  ERR  was  minimized  by  combining  the  Y(J)  values  for  all 

14 

vertical  lines. 


6.  SIMULATION  OF  CORRELATED  HORIZONTAL  COMPONENTS 


As  discussed  above,  the  correlations  in  the  horizontal  LOS  direction,  and  transverse 
direction,  t,  were  applied  to  the  previous  vertically  correlated  simulated  data  by  using  two- 
dimensional  Fourier  transform  synthesis.  The  model  autocorrelation  function  used  in  this 
report  leads  to  an  isotropic  two-dimensional  power  spectral  density  function  given  by: 


Haykin,  Simon,  ed.  (1991)  Advances  in  Spectnim  Analysis  and  Array  Processing,  VI, 
Prentice  Hall,  Englewood  Cliffs,  New  Jersey,  pp.  155-156. 
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2  2v 

a  va 


I  2  ,  2  ,  , 

nici  + 


for  -  oo  <  kpk  f<°°. 


This  representation  replaces  the  one-dimensional  -(2u+l)  power  dependence  with  a  two- 
dimensional  -(2u+2)  power  dependence.  was  used  to  filter  the  vertically  correlated 

Gaussian  random  values  at  a  given  altitude.  First,  the  procedure  found  the  two-dimensional 
Fast  Fourier  Transform  of  the  two-dimensional  Gaussian  correlated  values  at  a  given  altitude. 
These  spatial  frequency  values  then  were  multiplied  by  the  filter  function  -^jp(kf,k(J/AtA£.  A 

two-dimensional  inverse  Fourier  Transform  of  the  result  gave  the  final  correlated  simulated 
data  at  the  given  altitude.  The  process  stepped  up  in  altitude  to  yield  the  three-dimensional 
database.  Summarizing,  the  simulated  data,  S(t,  f),  for  a  layer  is  described  by. 

S(tJ)  =  FFT 

simulated  data  for  the  layer.  As  explained  in  the  appendix,  sampling  theory  results  caused  us  to 
modify  the  expressionfor  S(t,i). 


theoT 


AtAi 


X  FFT(G(ef,e^)).  where  G(et,e^)  are  the  vertically  correlated 


7.  RESULTS  -  VALIDATION  OF  SIMULATED  DATABASE 


A  practical  and  efficient  three-dimensional  non-stationary  stochastic  model  of 
atmospheric  structure  consistent  with  having  pre-assigned  vertical  and  horizontal  power 
spectral  densities  and  autocorrelation  functions  is  desired.  In  this  section,  we  examine  and 
compare  the  prescribed  structure  to  that  obtained  from  the  stochastic  hybrid  method  described 
earlier.  This  examination  serves  to  validate  the  database.  Correlation  lengths  and  o 
variances  were  taken  from  the  data  of  Strugala,  et.al.i^  As  discussed  above,  fitted  forms  of 
these  data  were  used  as  inputs  to  the  autoregressive  and  Fourier  transform  filters. 

First,  we  examine  the  power  spectral  densities  of  the  vertical  structure  as  depicted  in  the 
left  panel  of  Figure  6.  The  light  gray  solid  curve  depicts  a  one-dimensional  “theoretical” 
vertical  PSD  plotted  against  wavenumber  for  an  altitude  of  50.6  km.  This  figure  is  typical  of 
the  graphs  in  this  section  that  show  log-log  plots  of  PSDs  measured  in 

(8  temperature  I  temperature  f_  ^^venumber  measured  in  km-i.  At  50.6  km  the  input 

wavenumber 

parameters  have  a  vertical  correlation  length  Lcv=  1.72  km  and  a  vertical  spectral  slope 
(S)  =  -5/3.  The  constant  vertical  data  spacing  was  0. 1  km.  This  curve  may  be  compared  with 
the  PSD  of  the  AR  model  that  was  used  to  generate  the  vertically  simulated  data  and  with  the 
PSD  of  the  simulated  data  itself.  The  dashed  curve  in  this  plot  shows  the  one-dimensional 
vertical  PSD  of  the  autoregressive  model  that  was  used  at  the  50.6  km  altitude.  In  applying  the 
autoregressive  model,  the  variance,  o^,  was  normalized  to  a  value  of  1.  The  actual  variance 
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was  Introduced  into  the  synthesis  during  the  horizontal  computations.  A  6-coefficient  model 
was  used  to  calculate  or  “predict”  the  next  vertically  correlated  data  point  for  the  next  higher 
layer.  Points  in  the  layer  were  computed  by  applying  the  autoregressive  model  to  each  vertical 
line  in  the  array.  Shown  on  the  same  plot  is  the  PSD  computed  from  the  vertically  constructed 
data  around  50.6  km  altitude.  This  PSD  results  from  computing  a  6-coefficient  forward  and 
backward  autoregressive  model  of  the  vertically  simulated  data  from  layers  extending  from 
50.0  km  to  50.6  km  in  the  forward  direction  and  layers  extending  from  51.2  km  to  50.6  km  in 
the  backward  direction.  Comparison  of  the  curves  show  excellent  agreement  between  the 
prescribed  PSD  and  the  PSD  computed  from  the  synthetic  vertical  data. 

We  now  examine  the  horizontal  transverse  structure  as  depicted  in  the  center  panel  of 
Figure  6.  A  light  gray  solid  curve  depicts  a  one-dimensional  “theoretical”  transverse 
horizontal  PSD  plotted  against  wavenumber  for  an  altitude  of  50.6  km.  The  transverse 
direction  is  perpendicular  to  the  line-of-sight  direction  and  the  transverse  data  spacing  is 
approximately  0.100  km.  At  50.6  km  the  input  parameters  have  a  horizontal  correlation 
length  32.9  km,  and  the  variance  of  the  relative  temperature  fluctuations  =  8.34x10"'^. 
The  figure  represents  the  1-D  component  of  the  2-D  horizontal  PSD  having  a  high  frequency 
horizontal  spectral  slope  (S)  =  -8/3.  This  curve  may  be  compared  with  a  periodogram  of  the 
simulated  data  aligned  in  the  layer  along  the  transverse  direction.  It  may  also  be  compared 
with  a  PSD  computed  from  an  AR  model  of  that  same  data.  The  solid  curve  in  this  panel  shows 
a  one-dimensional  periodogram  of  the  transverse  simulated  data  in  the  layer  at  50.6  km.  The 
periodogram  has  a  Nyquist  frequency  of  approximately  l/(2x0. 1)  =  5  km‘i.  The  periodogram 
was  computed  by  averaging  the  individual  periodograms  of  the  simulated  data  along  each  of 
the  192  transverse  lines  in  the  layer.  Comparison  of  the  gray  and  solid  curves  shows  excellent 
agreement  between  the  prescribed  PSD  and  the  periodogram  computed  from  the  synthetic 
transverse  data.  The  dotted  curve  in  this  panel  shows  another  one-dimensional  transverse 
PSD  plotted  against  wavenumber  for  50.6  km  altitude.  The  curve  results  from  computing  a  10- 
coetficient  autoregressive  model  of  the  transverse  simulated  data  in  the  layer  at  50.6  km. 
Comparison  of  the  gray  and  dotted  curves  again  shows  excellent  agreement  between  the 
prescribed  PSD  and  the  PSD  computed  from  the  synthetic  transverse  data. 

Next,  we  examine  the  horizontal  line-of-sight  structure  as  depicted  in  the  right  panel  of 
Figure  6.  A  light  gray  solid  curve  depicts  a  one-dimensional  theoretically  aliased  line-of-sight 
(LOS)  horizontal  PSD  plotted  against  wavenumber  for  an  altitude  of  50.6  km.  For  a  detailed 
discussion  of  the  aliasing  issue  the  reader  is  referred  to  the  appendix.  The  LOS  data  spacing  is 
approximately  15.7  km.  The  curve  represents  the  1-D  component  of  the  2-D  horizontal  PSD 
having  a  horizontal  spectral  slope  (S)  =  -8/3.  This  figure  may  be  used  for  comparison  with  a 
periodogram  of  the  simulated  data  aligned  in  the  layer  along  the  line-of-sight  direction.  It 
may  also  be  used  for  comparison  with  a  PSD  computed  from  an  AR  model  of  that  same  data. 
The  solid  curve  in  this  panel  shows  a  one-dimensional  periodogram  of  the  line-of-sight 
simulated  data  in  the  layer  at  50.6  km.  The  periodogram  has  a  Nyquist  frequency  of  approxi¬ 
mately  l/(2xl5.7)  =  0.032  km‘^  The  periodogreun  was  computed  by  averaging  the  individual 
periodograms  of  the  simulated  data  along  each  of  the  512  LOS  lines  in  the  layer.  Comparison 
of  the  gray  and  solid  curves  shows  excellent  agreement  between  the  prescribed  PSD  and  the 
mean  of  the  periodogram  computed  from  the  synthetic  transverse  data.  The  reason  for  the 
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"noisy”  PSD  is  the  high  degree  of  correlation  between  neighboring  and  nearby  data  points. 

This  is  illustrated  in  Figure  7,  which  shows  the  PSDs  of  simulated  data  in  the  layer  along  the 
first  four  lines  along  the  LOS.  The  small  differences  between  the  PSDs  are  due  to  the  high 
degree  of  correlation  between  data  in  adjacent  lines.  The  dotted  curve  in  the  right  panel  of 
Figure  6  shows  another  one-dimensional  line-of-sight  PSD  plotted  against  wavenumber  for 
50.6  km  altitude.  This  curve  results  from  computing  a  10-coefficient  autoregressive  model  of 
the  line-of-sight  simulated  data  in  the  layer  at  50.6  km.  Comparison  of  gray  and  dotted  curves 
again  shows  excellent  agreement  between  the  prescribed  PSD  and  the  PSD  computed  from  the 
synthetic  line-of-sight  data. 

Power  spectral  density  computations  for  119.3  km  altitude  are  shown  in  Figure  8.  The 
graph  repeats  Figure  6  except  that  the  input  parameters  are  altitude  =  1 19.3  km,  =  0.0277, 

=  2.42  km,  and  =  49.1  km.  The  transverse  spacing  is  0. 101  km  and  LOS  spacing  is  15.9 
km.  Comparison  with  Figure  6  indicates  that  the  hybrid  method  described  herein  may  be 
employed  for  a  wide  range  of  altitudes.  Additionally,  -3  spectral  slope  models  have  been  run 
with  results  veiy  similar  to  the  -5/3  spectral  slope.  That  is.  the  power  law  dependence  and  PSD 
amplitudes  specified  by  the  models  are  preserved  in  the  simulated  data. 
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transverse  spacing 


Estimates  of  the  correlation  lengths  and  variances  for  the  individual  specifications  and 
realizations  presented  in  this  report  are  shown  in  Table  2.  Also,  Table  3  provides  the 
estimated  error  from  the  specified  parameter  values.  The  vertical  database  correlation  length 
estimates  are  derived  from  the  6-coefficient  forward  and  backward  estimates  of  the  vertical 
PSD  while  the  horizontal  database  correlation  length  estimates  are  derived  from  the  10- 
coefficient  estimates  of  the  horizontal  PSD  along  the  line-of-sight  direction.  The  estimated 
database  variances  are  the  actual  variances  computed  from  the  database  values  themselves  for 
the  specified  layer.  Examination  of  Tables  2  and  3  tends  to  indicate  that  the  error  of  an 
individual  realization  of  correlation  length  falls  within  about  30  percent.  Since  the  database 
represents  a  single  sample  from  an  ensemble  of  potential  synthetic  atmospheres  (different 
random  number  seeds  would  yield  additional  samples),  the  estimated  error  seems  reasonable. 
As  the  number  of  realizations  expand  it  is  expected  that  the  estimate  of  the  mean  would 
converge  to  the  specified  values. 


Table  2.  Correlation  Lengths:  Specified  vs  Database  Estimates 


Database  vs  Specified  Model  Parameters 


-5/3  slope 

-3  slope 

2 

o 

Lev 

Lch 

2 

o 

Lev 

Leh 

50  km 

spec. 

0.00083 

1.72 

32.9 

0.00083 

1.72 

32.9 

dbase 

0.00084 

1.53 

46.0 

0.00082 

2.17 

30.1 

89  km 

spec. 

0.0100 

1.80 

36.2 

0.0100 

1.80 

36.2 

dbase 

0.0092 

1.76 

33.1 

0.0096 

2.15 

29.0 

120  km 

spec. 

0.0277 

2.42 

49.1 

0.0277 

2.42 

49.1 

dbase 

0.0289 

1.62 

38.9 

0.0289 

2.21 

35.2 
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Table  3.  Estimated  Errors,  Specified  Parameters  vs  Database  Values 
Percent  Error,  Database  vs  Specified  Model  Parameters 


-5/3  slope 

-3  slope 

2 

a 

Lev 

Lch 

2 

a 

Lev 

Lch 

50  km 

1 

-11 

39 

H 

26 

-9 

89  km 

m 

-2 

-9 

■1 

19 

-20 

120  km 

m 

-33 

-21 

■1 

-9 

-28 

Figure  9  depicts  the  probability  distribution  histogram  of  the  simulated  data  in  the 
horizontal  layer  at  50.6  km,  while  Figure  10  depicts  the  cumulative  distribution  function. 
Visual  comparison  of  the  histogram  with  the  theoretical  (solid)  curves  indicate  that  the 
simulated  data  retains  the  Gaussian  shape  and  that  the  standard  deviation  of  the  data  falls 
within  8  percent  of  the  theoretical  value. 

8.  DATABASE 

A  random  access  data  file  consisting  of  68,91 1,104  data  points  (701  vertical  by  512 
horizontal  transverse  by  192  LOS)  was  constructed.  Each  point  was  saved  as  an  8  bit  byte, 
encoded  between  a  minimum  value  of  -128  and  a  maximum  value  of  +127.  A  small  separate 
corollary  ASCII  file  was  also  saved  that  provided  specification  information  about  the  random 
access  file.  The  ASCII  file  also  contained  sets  of  linear  transformation  constants  for  each 
altitude.  These  constants  allowed  efficient  computation  of  the  relative  temperature 

St 

fluctuations  from  the  byte  data  as  —  =  Ci  -vcn  xByteValue.  A  sample  of  the  ASCII  file  is 

T 

shown  in  Figure  1 1 . 


25 


UJ 


n 


ffi  o  «  8  ' 
CD  O  Ci  UJ  ‘ 

O  m  =3  CM  I 
2  Wzo( 
X  ^  lu  r 

g  t  g  ^1 

UJ  g  z  “  : 

— J  o  ^  . 

n:  2  X  <  : 


io-36z-t^  AiiSNsa  Ainigvaodd 


26 


Figure  9.  Histogram  Representing  the  Probability  Distribution  Function  of  the  Simulated  Data  at 
50.6  km.  Solid  line  is  the  theoretical  Gaussian  pdf.  The  variance  of  the  data  is  8.35x10''^. 


12345  =RAN0OM  NUMBER  SEED 
7-FEB-94  OUTPUT  FROM  PROGRAM  SCHMOD5 
5.0000E+01  THIS  IS  Z1  LOWEST  ALTITUDE  SIMULATED 
6.3700E+03  THIS  IS  RE  THE  RADIUS  OF  THE  EARTH 
TO  CALCULATE  THE  TRANSVERSE  DISTANCE  DT  AT  ALTITUDE  Z 
FOR  POINT  NUMBER  NPT,WITH  SPACING  AT  Z1  OF  DX 
(DIFFERENT  ALONG  THE  LINE  OF  SIGHT  AND  PERPENDICULAR 
TO  THE  LINE  OF  SIGHT), 

USE  THE  FORMULA  DT= ( RE+Z ) *DX* (NPT-1 )/( RE+Zl ) 

2.6667E+00  SLOPE  OF  PSD  HORIZONTAL  (2-D) 

1.6667E+00  SLOPE  OF  PSD  VERTICAL  (1-D) 

l.OOOOE-01  SPACING  (KM)  HORIZONTAL  RESULT  PERPENDICULAR  L.O.S.  LOWEST  ALTITUDE 
l.OOOOE-01  SPACING  (KM)  VERTICAL  RESULT 
512  NUMBER  OF  VALUES  PERPENDICULAR  LINE  OF  SIGHT 
701  NUMBER  OF  VERTICAL  VALUES 
192  NUMBER  OF  VALUES  ALONG  LINE  OF  SIGHT 
51  NUMBER  OF  SHEETS  FOR  INITIALIZATION 
1.5707E+01  SPACING  (KM)  ALONG  LINE  OF  SIGHT 
157  NUMBER  OF  CYCLES  FOR  ALIASING  L.O.S 

1  NUMBER  OF  CYCLES  FOR  ALIASING  PERPENDICULAR  L.O.S 
-5.7000E-01  =PSD  BIAS 
CREATE  A  FILE  CONSISTING  OF  BYTES 

EACH  OF  WHICH  REPRESENTS  A  NUMBER  BETWEEN  -128  AND  +127. 

THIS  IS  DIVIDED  INTO  RECORDS  AS  THE  FOLLOWING  LOGIC  INDICATES 
NUM=0 

DO  1=1, NUMBER  OF  ALTITUDES 

DO  J=l, NUMBER  OF  POINTS  ALONG  TO  THE  LINE  OF  SIGHT 
NUM=NUM+1 

WRITE  (8,REC=NUM)  (VALUE( I, K,J) ,K=1, NUMBER  OF  POINTS  PERPENDICULAR  L.O.S.) 

END  DO 
END  DO 

FOR  1=  1,ALTITUDE=  5 . OOOOE+OIKM,  (DT/T)=  -2.8273E-03+  7 . 8602E-04 *VALUE 

FOR  1=  2,ALTITUDE=  5 . OlOOE+OlKM,  (DT/T)=  -9.5905E-03+  7 . 8048E-04*VALUE 

FOR  1=  3,ALTITUDE=  5 . 0200E+01KM,  (DT/T)=  -1.1463E-02+  7 . 8550E-04 *VALUE 

FOR  1=  4,ALTITUDE=  5 . 0300E+01KM,  (DT/T)=  -1.2085E-02+  7 . 1519E-04*VALUE 

FOR  1=  5,ALTITUDE=  5 . 0400E+01KM,  (DT/T)=  -9.6289E-03+  7 . 2735E-04 *VALUE 

FOR  1=  6,ALTITUDE=  5 . 0500E+01KM,  (DT/T)=  -1.5482E-02+  8. 2105E-04*VALUE 

FOR  1=  7,ALTITUDE=  5 . 0600E+01KM,  (DT/T)=  -1.1941E-02+  8 . 7036E-04 *VALUE 

FOR  1=  8,ALTITUDE=  5 . 0700E+01KM,  (DT/T)=  -1.5864E-02+  8 . 8468E-04 *VALUE 

FOR  1=  9,ALTITUDE=  5 . 0800E+01KM,  (DT/T)=  -1.2681E-02+  9 . 0684E-04*VALUE 

FOR  1=  10,ALTITUDE=  5 . 0900E+01KM,  (DT/T)=  -8.8909E-03+  9 . 6629E-04*VALUE 

FOR  1=  11,ALTITUDE=  5 . lOOOE+OlKM,  (DT/T)=  -6.6515E-03+  8 . 9558E-04 *VALUE 

FOR  1=  12,ALTITUDE=  5 . IIOOE+OIKM,  (DT/T)=  3.0155E-03+  9 . 3763E-04*VALUE 

FOR  1=  13,ALTITUDE=  5 . 1200E+01KM,  (DT/T)=  1.3740E-04+  8 . 5331E-04 *VALUE 

FOR  1=  14,ALTITUDE=  5 . 1300E+01KM,  (DT/T)=  -1.4131E-03+  7 . 9499E-04 *VALUE 

FOR  1=  15,ALTITUDE=  5 . 1400E+01KM,  (DT/T)=  -9.1715E-03+  8 . 3941E-04*VALUE 

FOR  1=  16,ALTITUDE=  5 . 1500E+01KM,  (DT/T)=  -8.3961E-03+  9 . 5934E-04 *VALUE 

FOR  1=  17,ALTITUDE=  5.1600E+0im,  (DT/T)=  -1.4722E-03+  8 . 8723E-04*VALUE 

FOR  1=  18,ALTITUDE=  5 . 1700E+01KM,  (DT/T)=  -4.9128E-03+  8 . 2328E-04*VALUE 

FOR  1=  19,ALTITUDE=  5 . 1800E+01KM,  (DT/T)=  -4.3047E-03+  8 . 3197E-04 *VALUE 

FOR  1=  20,ALTITUDE=  5 . 1900E+01KM,  (DT/T)=  -6.3201E-03+  7 . 6369E-04*VALUE 

FOR  1=  21,ALTITUDE=  5 . 2000E+01KM,  (DT/T)=  2.5905E-03+  7 . 7723E-04 *VALUE 

FOR  1=  22,ALTITUDE=  5 . 2100E+01KM,  (DT/T)=  -4 . 8177E-03+  8 . 6719E-04*VALUE 

FOR  1=  23,ALTITUDE=  5 . 2200E+01KM,  (DT/T)=  -2.3910E-03+  9 . 4519E-04*VALUE 

FOR  1=  24,ALTITUDE=  5 . 2300E+01KM,  (DT/T)=  -1.0030E-02+  8 . 6416E-04 *VALUE 

FOR  1=  25,ALTITUDE=  5 . 2400E+01KM,  (DT/T)=  -1.4973E-02+  9 . 4909E-04*VALUE 

FOR  1=  26,ALTITUDE=  5 . 2500E+01KM,  (DT/T)=  5.9302E-04+  8 . 4785E-04 *VALUE 

FOR  1=  27,ALTITUDE=  5 . 2600E+01KM,  (DT/T)=  -4.3746E-03+  8 . 0043E-04 *VALUE 

FOR  1=  28,ALTITUDE=  5 . 2700E+01KM,  (DT/T)=  -1.9294E-03+  8 . 1456E-04*VALUE 

FOR  1=  29,ALTITUDE=  5 . 2800E+01KM,  (DT/T)=  -4.9838E-03+  8 . 1147E-04 *VALUE 


Figure  11.  Sample  Page  From  the  Database  Explanation  File. 
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9.  SYNTHETIC  2-D  SCENE  FOR  AN  IMAGINARY  TEMPERATURE  SENSOR 


The  availability  of  a  three-dimensional  temperature  structure  database  allows  us  to 
perform  computer  simulations  of  particular  sensors  viewing  such  an  atmospheric 
environment.  This  ability  enables  us  to  perform  economical  tests  and  experimentation  of 
more  fundamental  processes.  To  demonstrate  the  utility  of  the  database,  we  will  examine  the 
effects  of  viewing  the  simulated  atmosphere  from  an  imaginary  temperature  sensor  platform 
located  at  an  altitude  of  85  km.  We  will  use  a  database  generated  with  a  -5/3  power  spectral 
density  power  law.  For  these  purposes,  we  imagine  that  our  chosen  sensor  has  the  ability  to 
integrate  temperature  fluctuations  along  its  line-of-sight.  Our  computerized  arbitrary  sensor  is 
an  array  of  512  horizontal  pixels  by  350  vertical  pixels.  The  sensor  viewing  geometry  is  a 
rectangular  cone  having  a  full  horizontal  included  angle  of  2.19°  and  a  full  vertical  included 
angle  of  1.48°.  Thus  each  pixel  “looks”  along  its  own  line-of-sight  and  "sees”  0.0042°  in  both 
the  horizontal  and  vertical  direction.  The  sensor  plane  is  tilted  downward  with  the  center 
pixel  having  a  zenith  angle  of  95.23°.  This  geometry  prescribes  that  the  center  pixel  in  the 
bottom  row  “looks”  along  a  ray  that  has  a  tangent  altitude  of  50  km.  Since  we  wish  to  “view” 
as  much  of  the  simulated  atmosphere  as  possible,  the  forward  distance  at  85.2  km  altitude 
(away  from  the  sensor)  is  specified  to  have  an  arc  length  of  1,343  km.  At  this  distance,  the 
vertically  viewed  plane  perpendicular  to  the  line-of-sight  has  a  horizontal  and  vertical  spatial 
resolution  of  0.1  km.  The  geometry  is  shown  in  Figure  12. 

The  integrated  temperature  structure  along  selected  pixel  rows  in  zenith  look  angles  from 
the  minimum  of  94.6°  to  the  maximum  of  96°  is  depicted  in  Figure  13.  The  figure  shows 
computations  of  J{At(z)}dz  as  plotted  against  the  full  range  of  horizontal  look  angles  across 

the  sensor.  Here  At(z)  are  the  temperature  fluctuations  along  a  particular  line-of-sight  z. 
Figure  14  shows  similar  plots  but  this  time  the  temperature  structure  along  selected  pixel 
columns  comprising  horizontal  look  angles  from  -1.07°  to  +1.07°  is  plotted  against  the  full 
range  of  zenith  look  angles.  Taken  together,  the  plots  provide  a  two-dimensional  spatial 
visualization  of  the  integrated  horizontal  and  vertical  structure  as  “detected”  by  our 
imaginary  temperature  sensor.  These  spatial  series  may  be  analyzed  by  Fourier  Transform 
analysis  emd,  separately,  by  autoregressive  analysis  to  provide  estimates  of  the  power  spectral 
densities  of  the  integrated  structure  as  detected  by  the  sensor.  Figure  15  depicts  10  Fourier 
periodograms  of  the  horizontally  scanned  simulated  atmosphere  for  the  same  10  zenith  angles 
as  Figure  13.  Figure  16  uses  6-coefficient  autoregressive  analysis  to  provide  the  PSD  estimates. 
Examination  of  these  plots  indicate  that  the  power  spectral  density  slope  was  transformed 
from  the  -5/3  database  slope  to  an  image  plane  slope  of  -2.  This  is  a  surprising  and  not  well 
understood  result,  which  may  be  important  in  analyzing  atmospheric  radiance  measurements 
or  when  simulating  radiance  structure  from  correlation  overlay  maps.  Additional  two- 
dimensional  temperature  structure  scenes  were  computed  using  the  -3  database  spectral  slope. 
Results  showed  that  the  larger  slope  was  not  transformed:  instead,  the  image  plane  slope 
retained  the  -3  value.  Investigations  are  continuing  to  identify  the  processes  leading  to  these 
results. 
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ALTITUDE  OF  OBSERVER=  8.S0E+01  LOWEST  ALTITUDE  SIMULATED  =  S.OOE+01  2-D  SLOPE  OF  PSD  HORIZONTAL  =  2.67E  +00 

1-D  SLOPE  OF  PSD  VERTICAL=  1.67E+00  TRANSVERSE  SPACING  AT  LOWEST  ALTITUDE=  1,00E-01  VERTICAL  SPAC1NG=  1.00E-01 

SPACING  ALONG  LO.S.  =  1.57E+O1  NUMBER  TRANSVERSE  POINTS  =  512  NUMBER  LO.S,  P01NTS=192  NUMBER  OF  VERTICAL  SHEETS  =  701 

MAXIMUM  ANGLE  L.O.S.  USED  =  1.19E+01  DEGREES  NUMBER  OF  ZENITH  VALUES  CALCUUTED=350  CENTER  ZENITH  VALUE=  9.52E+01  DEGREES 

ZENITH  SPACING=  4.20E-03  DEGREES  NUMBER  OF  ALPHA  VALUES  CALCULATED=512  ALPHA  SPACING=  4.20E-03 
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Figure  13.  Integrated  Temperature  Structure  Along  Selected  Pixel  Rows  Comprising  Zenith  Look 
Angles  From  the  Minimum  of  94.6°  to  the  Maximum  of  96°. 
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Figure  14.  Integrated  Temperature  Structure  Along  Selected  Pixel  Colunms  Comprising  Horizontal  Look 
Angles  From  -1.07°  to  +1.07°. 
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Figure  15.  Fourier  Periodograms  of  the  Horizontally  Figure  16.  PSD  Estimates  of  the  Horizontally 

Scanned  Simulated  Atmosphere  for  the  Same  Look  Scanned  Simulated  Atmosphere  From  Auto- 

Angles  as  Figure  13.  regressive  Analysis. 


The  full  two-dimensional  integrated  temperature  structure  “scene”  is  depicted  in  Figure  17. 
Clearly,  the  figure  shows  that  the  three-dimensional  hybrid  method  provides  matrix 
realizations  having  significant  utility  for  computer  sensor  simulations.  The  technique 
preserves  relatively  large  horizontal  correlations  and  relatively  small  vertical  correlations. 
As  Figure  17  also  shows,  use  of  the  matrix  realizations  leads  to  a  better  understanding  of 
underlying  processes  and  provides  some  unexpected  results  in  the  two-dimensional  scene. 


10.  CONCLUSION 


Geophysical  phenomena  within  a  specified  domain  often  are  characterized  by  smooth 
continuous  power  spectral  densities  having  a  negative  power  law  slope.  A  stochastic  one¬ 
dimensional  autoregressive  approach  was  employed  to  generate  vertically  correlated  synthetic 
data,  and  two-dimensional  fast  Fourier  transform  synthesis  was  used  to  generate  horizontally 
correlated  data  from  the  vertically  structured  arrays.  The  joint  goals  of  reducing  the 
computational  burden  and  of  generating  three-dimensional  non-stationary  synthetic  structure 
databases  that  are  faithful  to  the  theoretical  descriptions  were  achieved.  A  complete  random 
access  database  file  consisting  of  68,911,104  points  was  generated  on  the  Phillips  Laboratory 
model  210  Convex  computer  with  an  average  execution  time  of  2.5  hrs.  The  process  preserved 
the  power  spectral  density  law,  correlation  scale,  variance,  and  probability  density  function.  A 
6-coefficient  autoregression  model  accurately  simulated  the  701  point  vertical  structure  and  a 
512  X  192  point  2-D  FFT  accurately  provided  the  horizontal  structure.  The  three-dimensional 
hybrid  method  provides  matrix  realizations  having  significant  utility  for  computer  sensor 
simulations.  We  have  demonstrated  the  utility  of  the  three-dimensional  database  and  we 
examined  the  effects  of  viewing  the  simulated  atmosphere  from  an  imaginary  temperature 
sensor  platform  located  at  em  altitude  of  85  km.  Use  of  the  matrix  realizations  leads  to  a 
better  understanding  of  underlying  processes  and  provides  some  unexpected  results  in  the  two- 
dimensional  scene.  We  have  shown  that  the  integrated  temperature  structure  preserves  the 
relatively  large  horizontal  correlations  and  relatively  small  vertical  correlations  but  that  the 
-5/3  power  spectral  density  slope  is  transformed  by  the  line-of-sight  integration  through 
several  layers.  This  is  a  surprising  and  not  well  understood  result  which  may  be  important  in 
analyzing  atmospheric  radiance  measurements  or  when  simulating  radiance  structure  from 
correlation  overlay  maps.  In  contrast,  the  -3  power  spectral  density  slope  was  unaltered  by  the 
integration  process.  Studies  are  continuing  with  the  hope  of  identifying  the  processes  leading 
to  these  different  results. 
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dimensional  Integrated  Temperature  Structure  "Scene' 
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Appendix 

Computational  Details 


Al. 


LINE-OF-SIGHT 


PSD  TRUNCATION  CORRECTION 


Since  the  constructed  database  has  a  finite  and  relatively  large  line-of-sight  data  spacing, 
a  natural  truncation  of  the  line-of-sight  power  spectral  density  function  must  occur  that  would 
cause  significant  aliasing  of  the  line-of-sight  PSDs  for  any  particular  realization.  The 
database  has  a  horizontal  transverse  Nyquist  frequency  about  5  km'^  and  a  line-of- 

sight  Nyquist  frequency  of  about  0.032  km'^.  Our  treatment  assumes  the  PSD  =  0 

beyond  .  Since  significant  energy  exists  between  j  and  .  aliasing  must  occur  in 

the  LOS  direction.  Consequently,  computational  “aliasing”  of  the  line-of-sight  PSDs  was  added 
to  the  “theoretical"  PSDs  by  the  following  formulation.  For  an  isotropic  two-dimensional 
power  density  spectrum  having  frequencies  k^  and  kg,  with  Nyquist  frequencies  k^^^  and 

{^NYg )  •  estimated  power  spectral  density  (PSDg)  is: 
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PSD^{kt,k^)  =  PSD  {kt,ki)  +  PSD  {kt,2kj^^-kf)  +  PSD  (k^ak^^  +k^) 


+  PSD  ^k^,4k^^  -k^j  +  PSD  ^k^,4k^^  +/c^j  +  --' 

until  the  second  argument  becomes  greater  than  •  Each  term  in  the  series  on  the  right 
forms  an  aliasing  branch  as  depicted  in  Figure  Al.  Two  special  cases  apply. 

At  k^  =0,  PSD^{kt,0)  =  PSD  (k^,0)  +  PSD  {k^,2k^^)  +  PSD  {lCt,4JcjvyJ+-  -  ■ 

Also,  at  Ic^  =  ,  PSDg  ,  ICfjY^  )-PSD  {kt,2kfnr^)+PSD  ikt,3kj^^)  +  PSD  (kt,5k^^)+ 

Hence,  the  computation  of  PSDg(/Cf,/c^)  was  the  “theoretical”  horizontal  PSD  that  was  used  in 
construction  of  the  simulated  database.  The  above  applies  since  the  two-dimensional  discrete 
Fourier  Transform,  F(/Ct,/c^),  of  a  function  /(t,f)  is: 

forO<lCt,/c^<n. 

j=0  1=0 

and  since, 

^-U{2K-k()  ^  ^-Ukt  ^  J  ^-U{2n*kf)  ^  ^-itki  ^  ^^-Uk^  J  _ 

This  implies  that  the  PSD(lCf,2jr-/c^)  contributes  to  the  PSD  of  fiUi)  at  frequency  kf.,ki  and 
the  PSD(/Cf,27i+lc^)  also  contributes  to  the  PSD  of 
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FILE  PCHM0D6  7-FEB-94  GREAT  CIRCLE  APPROXIMATION  EARTH  RADIUS  =  6.370E+03 
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Figure  Al.  Alias  Branches  of  Estimated  "Theoretical"  LOS  PSD. 


A2.  RECTANGULAR  WINDOW  CORRECTION  -  TRUNCATION  OF  AUTOCORRELATION 
FUNCTION 


The  horizontal  two-dimensional  discrete  Fourier  Transform  used  in  the  syntheses  of  the 
database  results  in  a  well-known  rectangular  window  sampling  effect.  Assuming  negligible 
truncation  of  the  autocorrelation  function  in  the  line-of-sight  direction,  truncation  in  the 
transverse  direction  leads  to  sidelobe  leakage  of  the  frequency  response  function  of  the 
digitized  rectangular  window.  In  the  transverse  direction,  the  effects  of  leakage  can  be 
computed  as  the  convolution  of  the  “theoretical”  fSD(/c^,/c^)  with  the  square  of  the  magnitude 

of  the  frequency  response  of  a  digital  sine  function.  The  magnitude  of  the  digital  sine  function 
resulting  from  sampling  N  values  equally  spaced  by  At  in  the  transverse  direction  is: 


H(k,) 


At 


sin  (nkf^NAt) 
sin(jt/CfAt) 


This  means  the  PSD,  after  accounting  for  leakage  in  the  transverse  direction,  is  the  digital 
adjusted  function,  PSDgi(kf,fc^)  where. 


NAt 


•  2r 
sin  1 

iTtNAt 

.  2 
Sin 

[(kf  -/c)TtAt 

and  where  we  assume  PSDgi^|/Cf|>/c^^,/c^j  =  0.  Normalizing  so  there  is  no  change  in  the  total 

energy,  that  is  to  say  o^,  we  arrive  at  the  “corrected”  power  spectral  density  function, 

PSD*  The  normalization  is  necessary  because  the  variance  of  the  “continuous”  power 

spectral  density  differs  from  the  variance  of  the  “discretized”  power  spectral  density.  The 
normalization  depends  on  the  continuous  case  variance  ^  defined  by 

af(  =  f  PSD(k,k()dlc,  and  also  on  the  discrete  case  variance  ^  defined  by 


'■NYf 


1 


NAt 


1 


y  (with  step  size 

NAt  ,  ^  t  (}  F 

itct=-kNYf 
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The  "corrected”  PSD  then  becomes: 


PSD*  (kt  ,k^}  =  X  PSD®  (/Cf  ,k^). 


Accounting  for  the  discretized  windowing  effect,  the  two-dimensional  white  noise  spectrum 
should  be  multiplied  by  some  value  A  (kf,k^).  Assuming  negligible  contribution  from  the 
truncation  of  the  autocorrelation  function  in  the  line-of-sight  direction,  A  {k^,k()  is 
determined  by: 


A*{kt,k^): 


PSD*  {kt,ki) 


where  At  is  the  transverse  spacing  and  At  is  the  LOS  spacing.  Using  the  expression  for 
A* {kf.,k(),  a  computational  "correction”  to  each  LOS  frequency  may  be  made  to  account  for 
the  windowing  effect  but  such  computation  is  prohibitively  intensive.  Instead,  the  following 
approximate  scheme  was  devised. 

Since  by  construction,  f  |H(/c)|  dk  =  l.  the  convolution  theorem  requires 

NAt 

^2  n  * 

X  =  Observation  of  A  (k^.k^)  indicates  that  the  major  sampling  effect 

fcj  kf 

occurs  for  the  PSD(k^,k(}  at  the  extrema  =  0.  This  observation  is  illustrated  in  Figure  A2 
which  shows  a  "theoretical”  transverse  PSD  with  k^  =  0  versus  a  512  point  PSD  “corrected  ,  as 
described  above,  for  truncation  of  the  autocorrelation  function.  The  plot  has  the  same  input 
parameters  previously  given  for  50.6  km  altitude.  The  one-dimensional  high  frequency 
spectral  slope  is  (S)  =  -5/3,  =  8.34  x  lO'^,  L^^  =  32.9  km,  and  the  transverse  spacing  =  0. 1 

km.  The  plot  indicates  that  the  “corrected”  PSD  closely  follows  the  “theoretical”  PSD  for  all 
the  discrete  frequencies  except  the  zero  frequency.  To  avoid  the  time-consuming  convolution 
integral,  an  approximation  was  developed  that  modified  the  PSD(/Cf,lc^)  at  k^  -0.  The 


approximation  was  devised  from  the  following  reasoning.  We  wish  the  variance  of  our 
simulation  to  equal  a  given  value  but  the  actual  value  of  the  “uncorrected”  simulated 


variance  is  equal  to  ^  ^A^(k(.,k^j,  where  and  where  PSD^{kf^,k^)  is  the 


two-dimensional  theoretical  model  with  aliasing  in  the  k^  direction: 
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PSDg  (/Cf./c^)-  PSD(/Cf,/c^]  +  PSD  I  /Cf,2/Cjyy^  -kg'j  +  PSDlki.,2k^^  +/c^j 


+  PSD^fcj,  4fcjyy^  —k(j  +  PSD^kf^,4k^  ^  j  +  ■  ■  • . 

2 

Thus,  we  wish  to  find  a  modified  approximate  value,  A ,  such  that  the  theoretical  variance,  , 
is:  =  ^  ^A^(lCf,lc^).  As  noted,  with  the  extrema  at  fcj  =0,  we  set  a(0,/c^)  by  the  linear 

kt  kf 

interpolation  relation:  A(0,/c^)  =  gA(0,/c^)  +  (l-fif)A(lC(_2,/c^),  where  kt_-i  is  the  first  non-zero 
LOS  frequency  value  and  we  set  all  other  A(kf,k()  =  A(kf.,k().  Then  it  follows: 

kf^O 


k(  kf^Q  k( 


Expanding,  we  find  a  quadratic  expression  for  g,  which  is: 


0  =  g''J^(A^{0.k,)  +  A\kt=„k,)-2A{0,k,)A{kt^„k,)) 


+  9  X  (2  A  (0.  t, )  A  ()c„„  )c, )  -  2  a"  (lc,,„  Ic,)) 

kf 


kf  kf^Q  kf 


Figure  A3  illustrates  the  computation  of  g  on  a  log-log  plot  of  g  versus 

i'rt 

near  power  law  relationship  between  g  and  — 

At 


.  The  plot  shows  a 
At 
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T ransverse  PSD 


10-6  10-5  10-4  10-3  10-2  10-1  10°  10^ 
Wavenumber  (km-"') 


Figure  A2.  "Theoretical"  Transverse  PSD  Versus  a  512  Point  PSD  "Corrected"  for 
Truncation  of  the  Autocorrelation  Function.  The  plot  has  the  same  input  para¬ 
meters  previously  given  for  50.6  km  altitude.  The  one-dimensional  high  frequency 
spectral  slope  is  (S)  =  -5/3,  =  8.34  x  IQ-^.  =  32.9  km  and.  transverse 

spacing  =  0.1  km. 
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2-D  PSD  SLOPE=2.67E+00  SPACING  PERPEI^piCULAR  TO  L.O.S.  LOWEST  ALTITUDE-  1.00E-01  SPACING  ALONG  L.O.S.  LOWEST  ALTITUDE-  1.57E+01 
STRAIGHT  LINE  FIT—1.84E-01  +(LOG10(LCH)-  2.57E+00)* -3.42E-01  NUMBER  OF  ALIASING  VALUES-1 57 


Figure  A3.  Approximation  Correction  for  Rectangular  Window  Discrete  Fourier  Transform  Digital 
Sine  Function. 


